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ABSTRACT 

A parameterized model of the mass distribution within the Milky Way is fitted to the 
available observational constraints. The most important single parameter is the ratio 
of the scale length i?^,* of the stellar disk to Rq. The disk and bulge dominate v c (R) at 
R ^ Ro only for Rd,*/Ro ^ 0.3. Since the only knowledge we have of the halo derives 
from studies like the present one, we allow it to contribute to the density at all radii. 
When allowed this freedom, however, the halo causes changes in assumptions relating 
to R <C Rq to affect profoundly the structure of the best-fitting model at R ^> Rq. 
For example, changing the disk slightly from an exponential surface-density profile 
significantly changes the form of v c (R) at R ^> Rq, where the disk makes a negligible 
contribution to v c . Moreover, minor changes in the constraints can cause the halo 
to develop a deep hole at its centre that is not physically plausible. These problems 
call into question the proposition that flat rotation curves arise because galaxies have 
physically distinct halos rather than outwards-increasing mass-to-light ratios. 

The mass distribution of the Galaxy and the relative importance of its various 
components will remain very uncertain until more observational data can be used to 
constrain mass models. Data that constrain the Galactic force field at z *J R and at 
R > Rq are especially important. 
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1 INTRODUCTION 



. ^ One of the fundamental tasks of Galactic astronomy is the 
■ determination of the mass and luminosity distributions of 
S_t ' the Milky Way. In the 1950s the development of radio as- 
?3 t tronomy opened up the study of the Galaxy's large-scale 
structure, and much of the understanding that was attained 
at that time was summarized by Schmidt's (1956) mass mo- 
del. In the 1970s and early 1980s our picture of the Milky 
Way changed in response both to studies of external galaxies 
and to a growing awareness of the existence of "dark matter" 
at large radii. These developments were reflected in the Bah- 
call & Soneira (1980), Caldwell & Ostriker (1981) and Rohlfs 
& Kreitschmann (1988) Galaxy models. These models were 
based on two rather different methodologies: whereas Bali- 
call & Soneira concentrated on fitting the distribution of 
luminosity within the Galaxy by fitting star counts, Cald- 
well & Ostriker and Rohlfs & Kreitschmann concentrated on 
fitting various measures of the Galactic gravitational force- 
field. Never the less, all these models decompose the Galaxy 
into "components" that are motivated by photometric stud- 
ies of external galaxies, and incorporated a range of dynami- 
cal constraints. It is our aim in this and subsequent papers 
to update and extend these models. 

The principal direction in which we wish to extend tra- 
ditional galaxy models is the incorporation of kinematic in- 



formation that is capable of constraining the degree of flat- 
tening of the mass distribution. The kinematic information 
that has traditionally been used to constrain galaxy models 
- the shape of the circular-speed curve, the values of the Oort 
constants, etc - relates almost exclusively to the radial force 
within the plane. Such information is in principle incapable 
of determining how much of the Galaxy's mass lies near the 
plane, which is clearly of prime importance astrophysically. 

This deficit of kinematic information has been papered 
over in two ways. The first is a one-dimensional analysis of 
the vertical structure of the disk along the lines pioneered by 
Oort (1932). Such analyses ultimately come up against the 
problem that the vertical and horizontal motions of stars 
do not decouple to the necessary degree, so that a one- 
dimensional analysis cannot precisely determine the vertical 
distribution of matter - see, e.g., Kuijken & Gilmore (1991). 

The second way in which the deficit of kinematic con- 
straints on the Galaxy's vertical structure has been papered 
over is to use star counts to constrain the flattening of the 
components into which the overall luminosity density is de- 
composed, and to assume that each component i is charac- 
terized by a constant mass-to- light ratio Yj. This procedure 
is methodologically questionable since it appears to presume 
that the phenomenon of "dark matter" implies the existence 
of a completely dark component comprised of exotic parti- 
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cles, whereas it may merely reflect the variation from point 
to point of each TV 

This is the first of a series of papers in which we plan to 
overcome these difficulties by treating the orbits of stars in 
the meridional plane with sufficient sophistication. Our ap- 
proach, which has been described elsewhere (Binney f994, 
Dehnen & Binney f996), is iterative: we choose a potential, 
determine a range of orbits in this potential, populate these 
with stars of various spectral types and then compare the re- 
sulting predictions with the available surveys. The potential 
is modified in the light of this comparison. Thus our first step 
is to choose potentials (and thus, implicitly, mass models) 
which are compatible with all the standard kinematic and 
photometric constraints. This task of updating the Caldwell 
& Ostriker and Rohlfs & Kreitschmann models is even now 
by no means trivial, as is made apparent by recent diver- 
gent results (Gates, Gyuk & Turner, f995; Cowsik, Ratnam 
& Bhattacharjee, f996; Evans 1996). Therefore we believe 
it will be useful to present in this paper our initial mass 
models and the consider actions upon which they are based. 
Computer programs for evaluating the density and poten- 
tial of the models are available upon request to the authors. 
Although there is now abundant evidence that the inner Ga- 
laxy is significantly non-axisymmetric, our model conforms 
to the traditional axisymmetric pattern because the success 
of axisymmetric models in accounting for observations in 
the 21-cm line of hydrogen at longitudes I J> 30 deg sug- 
gests that orbits in the Galactic potential that carry stars 
to radii r <; 5kpc can be accurately modelled by orbits in 
an axisymmetric potential. 

We do not distinguish between the visible halo, of which 
RR-Lyrae stars and metal-poor globular clusters are clas- 
sical tracers, and the putative dark halo: since we do not 
understand why the mass-to-light ratio rises with galacto- 
centric radius r, we are at liberty to assume that the Galaxy 
possesses a single, massive halo that simply becomes more 
luminous with decreasing r. 



2 FUNCTIONAL FORMS 

Our mass model contains three principal components: the 
disk, the bulge and the halo. 



2.1 The Disk 

Our disk is made up of three components, namely the ISM, 
and the thin and thick stellar disks. The density of each 
sub- disk is given by 



, n \ E d / Rm R \z\ 
Pd(R, Z) = 7T- ex P W "p 

2z d \ R Rd z A 



(1) 



With Rm = 0, equation (|l|) describes a standard double 
exponential disk with scale- length Rd, scale-height Zd and 
central surface-density Ed. Since there appears to be very 
little interstellar gas between the molecular ring at R — 4- 
5kpc and the nuclear disk at < 200 pc (Dame et al. 1987), 
there should be a depression in the central surface-density of 
the ISM. The parameter R m in equation (|l|) allows for such 
a central depression. We set R m — for the stellar disks 
and adopt R m = 4 kpc for the ISM. The total mass of a disk 
with density ([!]) is 



Table 1. Fixed parameters of the disk components 

Component Contribution Rd/Rd.* Rm Zd 
to S(_R ) 



ISM 0.25 
thin disk 0.70 
thick disk 0.05 



2 4 kpc 40 pc 
1 180 pc 

1 1000 pc 



M d = 4tt E d Rm Rd Ki (2 V / R m /Rd) 



(2) 



where Ki is a modified Bessel function. For R m = this 
gives M d = 27rE d 7?d- 

Table |l| gives our adopted values of the scale heights 
of the three sub-disks as well as the fraction of the whole 
disk's surface density at Rq which is contributed by each 
sub-disk. As can also be seen from this table, we fix the ratios 
between the sub-disk's scale- lengths. The value of Rd., the 
scale length of the stellar disk, and the mass of the whole 
disk are obtained from least-squared fits to the observational 
constraints to be discussed below. 



2.2 Bulge and Halo 

The bulge and halo are each described by the spheroidal 
density distribution 

_2/_a 



(3) 



where 



(4) 



„, — I n>2 i -2 2 

m = (R + q z 

Thus the density of bulge and halo is proportional to r~ 7 
for r <C ro, proportional to for m < r < r c , and is 
softly truncated at r = r±. Infrared photometry obtained by 
the COBE/DIRBE satellite and analyzed by Spergel, Mal- 
hotra & Blitz (1997) yields values for four of the five bulge 
parameters: we adopt f3b — 7b = 1.8, <?b = 0.6, ro,b = 1 kpc, 
and r t ,b = 1.9 kpc. The density normalization po,b, which 
is not determined by the COBE/DIRBE data, is obtained 
from our least-squared fits. 

The axis ratio of the halo is not significantly constrained 
by the observations discussed below, and we arbritraily set 
it to t/h = 0.8 (our standard value). However, we will also 
consider models with = 0.3 motivated by (i) possible ev- 
i dence that t he halos of external galaxies might be that flat 
( |011ing 1996] Sackett et al. 1994), and (ii) Sciama's (1990) 



decaying neutrino model for dark matter that requires a 
high local neutrino density (~ 0.04Mqpc -2 ) only achiev- 
able with a flat neutrino halo. 

The remaining five halo parameters, /3h, 7h, po,h, r o,h, 
and rt,h, are determined by the least-squares fitting proce- 
dure described below subject to the restrictions — 2 < 7h < 
7b and /3h > 1, which limit the sharpness of any inner and 
outer edges of the halo, respectively. 



2.3 Gravitational potential and forces 

The total gravitational potential $ of a model must satisfy 
Poisson's equation 

,2 3 

V $ \ - \ - 

2_^P*,i + / ,P<M- 

j=l i=l 



4ttG 



(•>) 
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Table 2. Pairs of functions h(z) and H (z) defined in section 



h(z) = d 2 H/dz 2 


H{z) 


5{z) 


M 
2 


1 exp(-M) 


^[exp(-M)-l + M] 


-r— sech 2 7t— 
4z d 2z d 


In cosh t^— 



A standard way to solve (^) involves expanding p in spherical 
harmonics (c.f. BT section 2.8). Unfortunately, the expan- 
sion of a thin disk converges very slowly, so this straightfor- 
ward approach does not yield a fast and accurate solution 
of Poisson's equation. 

Fortunately, Kuijken & Dubinski (1994) have described 
a modified multipole technique that works well when the 
density is sum over components that are separable in cylin- 
drical coordinates, that is, are of the form 

p d ,i = fi(R)hi( Z ), (6) 

where 1 = f_ htdz so that /(ii) is the radial surface- 
density profile. Let Hi(z) be such that 

(i) Hi'(z) = h i (z) 
and (ii) Hi(Q) = Hj(0) = 0, 
where the prime denotes a derivative as usual. Then we write 



(7) 



(8) 



where the argument of fi is now the spherical radius r rather 
than the cylindrical radius R and <&me is a function to be 
determined. At z = the second term on the right hand side 
of (|8|) and its first derivatives vanish, so both the potential 
and the forces in the plane are determined by <&me alone. 
Inserting (|^) and (^) into (^) we obtain for <&me 



AtiG 



fl'{r)Hi{z) - -f^r) [H t {z) + zHi{z)] 



(9) 



This equation takes the form of Poisson's equation for $mb 
with a mass-density given by the complex expression on its 
right-hand side. At z = we have ii = r, so with ([?]) this 
expression simplifies to Ps,j- That is, $me is generated 
by a mass distribution that is not strongly confined to the 
plane, and can be economically evaluated by expanding both 
sides of equation (^) in spherical harmonics. 

Once $mb has been found, it and its first derivatives are 
stored on a grid in lnr and \z/r\. A two-dimensional fifth- 
order spline is used to interpolate on this grid: at each grid 
point this spline yields the stored values of the potential and 
its derivatives, and the forces have everywhere continuous 
first and second derivatives. In particular, the interpolated 
forces agree with the derivatives of the interpolated poten- 
tial, as is necessary if energy is to accurately conserved along 
numerically integrated orbits. Furthermore, the evaluation 
of potential and forces is quick once the spline coefficients 
have been computed. 

We plan to make available after publication of this pa- 
per a C++ source code for the evaluation of the potential 



of any superposition of spheroids (eq. pf) and exponential 
disks with vertical profiles as in Table g - send e-mail to 
w.dehnen@physics.ox.ac.uk. (Public-domain C compilers are 
available that allow C++ code to be linked to otherwise pure 
C programs.) 



3 THE OBSERVATIONAL CONSTRAINTS 

Three groups of observational data constrain the values of 
the free parameters in the model introduced in section ^| 
(i) tangent velocities at R < Ro; (ii) rotation velocities at 
R > Ro; (iii) other data, such as the values of Oort's con- 
stants and the local surface density. We discuss each group 
of constraints separately. 



3.1 Terminal velocities for the inner Galaxy 

For an axisymmetric Galaxy with circularly rotating ISM, 
the peak velocity along a given line-of-sight at b = and I in 
either the first or fourth quadrant originates from the radius 
R = Rosinl. Relative to the LSR this 'terminal velocity' is 
related to the circular speed v c by 



I'tc 



= v c (Ro sin/) — v c (Ro) sinL 



(10) 



In reality, non-circular motions of the ISM induced, for in- 
stance, by spiral arms, lead to deviations from this ideal 
relation. However, outside the region |Z| ^ 20deg that is 
dominated by the bar, these deviations are expected to be 
ignorable for our purposes. 

Numerous surveys of the ISM have been undertaken. 
I n this study we restrict ourselves to three surveys in H I 



flWeaver fc Williams 1973] IfBania fc Lockman 19841 Ker r 
et al. 1986) and one in CO (Knapp, Stark fc Wilson 1985|). 



Malhotra (1994,1995) has modelled these raw data in detail 
and kindly provided us her values for the terminal veloci- 
ties in electronically readable form. We restrict the data to 
I sin Z| > 0.3 to avoid distortions by the central bar. 

3.2 The rotation curve of the outer Galaxy 

For an axisymmetric galaxy, the radial velocity relative to 
the LSR, t?isr, of a circularly orbiting object at galactic co- 
ordinates (I, b) and galactocentric radius R is related to the 
circular speed by 



W(R) = . ^ = 2}.v e (R) - v c (R ) 
sin I cos R 



(11) 



As is well known, for R > Ro one cannot infer ii for an 
object at given / without a knowledge of the distance d to 
the object. If d is known, then ii follows from 

R = (d 2 cos 2 6 + R 2 ~ 2R d cos b cos l) 1/2 . (12) 

Several studies are available that contain measured va- 
lues of W and d for objects that ought to be on nearly 
circular orbits. Here we use the data of two recent stud- 
ies. Brand & Blitz (1993) list Hn regions/reflection nebulae 
that have (spectro-) photometric distances and associated 
molecular clouds with measured radial velocities. Pont et al. 
(1997) give radial velocities and photometry for classical 
cepheids in the outer disk. For these objects we have trans- 
formed v r to ui S r and evaluated the distances using the 
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period- luminosity relation derived from Hipparcos paralaxes 



by Fea st and co-wo rkers (Feast 
Whitel( fck~1997a,b| ), 



Catchpole 1997 



Feast & 



M v = -2.81 log P - 1.43, in conjunc- 
tion with a Galactic period-colour relation from Laney & 
Stobie (1994): (B - V)o = 0.416 log P + 0.314. 

We have rejected objects for which either 155° < / < 
205° , or d < 1 kpc, or W > 0, because for these objects v\ BT is 
very likely dominated by non-circular motions. Furthermore, 
we have not used data points at R < Ro where the terminal 
velocities provide a much better constraint. 93 of the 205 
objects in Brand & Blitz (1993) and 26 of the 48 Cepheids 
survived this cull. 

We have not employed a number of data sets that are 
similar in scope to those of Brand & Blitz and Feast et al. 
because these sets are either rather restricted in their ra- 
dial coverage, or have problematic distances. For example, 
the distances to carbon stars are seriously affected by both 
Malmquist bias and interstellar extinction (Schechter, pri- 
vate communication). 

A technique for measuring W(R) from 21-cm emission 
without independent distance information has been pro- 
posed by Merrifield (1992). This involves a determination 
of the extent in b of the emission observed at given W: in 
an axisymmetric galaxy with circularly orbiting H I, all emis- 
sion at given W will originate in a galactocentric ring. If this 
ring has a constant vertical extension, it creates a distinct 
pattern in the (/, b) plane. From the H I surveys of Weaver & 
Williams (1974) and Kerr et al. (1986) Merrifield estimated 
relative galactocentric distances R/Ro by fitting this char- 
acteristic pattern to the emission from each bin in W. It 
is not immediately apparent how accurate Merrifield 's radii 
are since both random motions in the plane and systematic 
variations in the thickness of the H I layer around circles will 
contribute errors. 

It is worth comparing the three data sets before try- 
ing to fit them to a mass model. Fig. [I] shows the data 
in W vs. R/Ro for Ro = 7.5 kpc, 8 kpc, and 8.5 kpc. Also 
shown is W(R) predicted by flat rotation curves for v c (Ro) = 
180 km s -1 (uppermost line), 200 km s -1 , 220 km s -1 , and 
240 km s -1 (lowest line). Several points can be drawn im- 
mediately from this Figure, (i) The data of Merrifield (us- 
ing Hi) at R > 1.3i?o are incompatible with the Cepheid 
data (Pont et al. 1997). (ii) The Hn data at R < 1.5R 
and the Cepheid data imply a falling rotation curve, un- 
less v c (Ro) ^ 240 kms -1 and/or R & 7.5 kpc. (iii) Merri- 
field's measurements require a rising rotation curve, unless 
Vc{Ro) £ 180 kms -1 . 

Clearly, it makes no sense trying to fit inconsistent data; 
we must decide which data to trust. Cepheid distances have 
been extensively studied, whereas Merrifield's novel method 
rests on several ad hoc assumptions. Furthermore, his values 
for W are suspect because they are non-monotonic at R > 
2Ro- (Note that they do not represent individual objects as 
do the other data points.) Therefore, we have decided to 
discard Merrifield's points altogether. 



3.3 Other constraints 



By dividing equation ( |10D by Ro sin I and equation (|llj) by 
Ro, one sees that studies of the ISM measure SI = v c (R)/R 
at various positions in the Milky Way relative to its local 




Figure 1. Data on the outer rotation curve: W vs. R/Ro (error- 
bars are omitted for clarity). Squares, stars, and triangles repre- 
sent Hn regions, Cepheids, and Hi measurements, respectively. 
The lines refer to flat rotation curves with v c (Rq) = 180kms — 1 
(solid), 200 km s" 1 (dotted), 220 km s" 1 (short dashed), and 
240 kms -1 (long dashed), respectively. Each panel corresponds to 
a different assumed value of Ro; increasing Ro shifts the squares 
and the stars to the left. 



value. To fix the absolute values of v c , additional information 
is essential. 



3.3.1 Oort's constants 

Oort's constants A and B are defined by 



A = 



B = 



1 /Vc _ OVc \ 

2 \R dRJ 



1 /Vc dv c \ 

2 \R dRJ 



(13) 



They can be derived from the kinematics of nearby stars. 
Kerr & Lynden-Bell (1986) reviewed the published measu- 
rements and concluded that A = 14.4±1.2, B = -12.0±2.8, 
and A — B — 26.4 ± 1.9, all in units of kms -1 kpc -1 . Han- 
son (1987) found from an analysis of the NPM catalog A = 
11. 3 ±1.1, B = -13.9 ±0.9, and A-B = 25.2 ±1.6, while a 
more recent study of Hipparcos proper motions for Cepheids 
by Feast & Whitelock (1997b) yields A = 14.8 ± 0.8, and 
A — B = 27.2 ± 0.9 in the same units. Since the Hipparcos 
values are likely to be significantly more reliable than the 
earlier studies, we give them the highest weight. 

A potentially relevant measurement is of the proper mo- 
tion of SgrA*: m = (-6.55 ± 17) masyr' 1 = (-31.05 ± 
0.81) kms -1 kpc -1 (Backer 1996). To this the peculiar mo- 
tion of the Sun contributes Vq/Rq — 1 kms -1 kpc -1 . Hence, 
if we believe tha t Sgr A*, which is certainly massive (Eckart 
&i Genzel 1997), is stationary at the centre of the Galaxy 
then we have A — B = 30 ± lkms -1 kpc -1 , which is sig- 
nificantly inconsistent with the result of Feast & Whitelock. 
In view of this conflict and the possibility that the Galactic 
centre is oscillating, perhaps in an m = 1 mode, (Gould & 
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Ramirez 1997) we have not used Backer's result but imposed 
the constraints 1 



A 



A= 14.5 ± 1.5 kms-'kpc" 1 
B = -12.5±2 kms-'kpc" 1 
B= 27 ± 1.5 kms-'kpc" 1 



(14) 



3.3.2 The mass at large radii 

While observations for A, B, and v c (Rq) restrict the circu- 
lar speed locally, and hence the mass inside Rq, there are 
some important constraints at much larger radii. The total 
mass within a sphere of radius r S> -Ro can be estimated 
(i) from the velocity distribution of the Milky Way's satel- 
lites, (ii) from the maximal locally observed stellar velocity 
(the 'escape velocity' argument), (iii) the timing of the local 
group, and (iv) by modelling the dynamics of the Magel- 
lanic Clouds and Stream. All these estimates rely on certain 
assumptions and are model dependent. However, with rea- 
sonable assumptions Kochanek (1996) found a simple model 
that satisfies the first three of these constraints and yields an 
acceptable value of v c (Ro). From his Figure 7 we extracted 
for the mass inside 100 kpc -Mfl<iookpc = 7.5 ±2 x 10 11 M Q . 
For comparison, by modelling the dynamics of the Magel- 
lanic Clouds and Stream Lin, Jones & Klemola (1995) found 
that Mfl<iook P c = 5.5 ± 1 x 10 11 Mq. Clearly, the uncertain- 
ties here are dominated by systematic errors, so we allow for 
a generous error bound on our adopted constraint, which is 

= (7 ±2.5) x 1O U M . (15) 



3. 3. 3 The local vertical force 

The vertical force K z at some height above the plane places a 
condition on the local mass distribution, and certainly is an 
important observable our model must agree with. Using K 
stars as a tracer population, Kuijken & Gilmore (1989,1991) 
have deduced 



\K x (Bo,l.l kpc) | =2ttGx (71±6)M oP c~ 2 . (16) 



We have adopted this as a constraint for our models. 



3.3.4 The disk's local surface density 

Unfortunately, the local disk surface density £0 is no t as well 
determined as th e closely related quantity K Zll 1 (Kuijken 



&c Gilmore 1991 ). However, by counting identified matter 
Kuijken & Gilmore (1989) concluded that £ s tars+gas(-Ro) = 
48 ± 8 Mq pc -2 . We adopt the constraint 



£0 = EdiskORo) >40 M Q pc- 



3.3.5 The dispersion velocity in Baade's window 



(17) 



Finally, an important constraint on the bulge is provided 
by the observed velocity di spersion of the bulge in Baade's 
window of 117±5kms _1 (Rich 1988; Terndrup, Sadler & 



|^| The product ARo is not an independent observable, as it is 
usually measured from the ISM's terminal velocities, which have 
already been employed as constraints. 



Rich 1995). The simplest way to estimate the velocity disper- 
sion of our models is to solve the Jeans equations assuming 
isotropy in the velocities, which yields 



2 



1 

Pb 



Pb -5- dz, 

az 



(18) 



where the subscript b stands for bulge. The bulge is now 
known to be significantly elongated towards us (e.g., Bin- 
ney, Gerhard & Spergel, 1997) and this elongation is prob- 
ably reflected in the line-of-sight dispersion along the Ga- 
laxy's minor axis being larger than equation dli| ) allows. On 
the other hand, the velocity dispersion probably falls as one 
moves away from the minor axis along the line of sight, and 
this effect will tend to cause equation (^) to overestimate 
the measured dispersion within Baade's window. In view of 
these oppositely directed factors, we adopt ([ll]) as the cen- 
tral value of our constraint on the dispersion within Baade's 
window and allow a wide range around this value: 

a BW = (j b (0.0175R ,-0.068i?o) = 117 ± 15kms -1 . (19) 



4 FITTING THE MASS MODEL 

The free parameters of the mass model described in section 
^| are determined by minimizing the quantity 



2 _ rv in 2 . "nut 2 
Xtot — jy. X in + M Xaut + 



2 

"^othe 



(20) 



that is the sum of pseudo-chi-squared contributions from 
our three classes of constraint. Here the Ni are the num- 
bers of data points actually used, while the Wi are weights, 
which may be interpreted as the number of really indepen- 
dent constraints (for instance, Oort's A has been obtained 
from much more data than we use for Vterm)- Clearly, the 
Wi are subject to ones prejudices, we took W in = W out = 
W otllet =N othet =6. 

There are contributions to xf n from 53 data points at 
I < and 77 at I > 0. In order to minimize the influence of 
systematic deviations from circular motion, which differ on 
the two sides of the Galaxy, the data for positive and nega- 
tive longitude are weighted by 0.844 and 1.23, respectively. 
This gives an effective number of 65 data points on either 
side, while leaving the effective total number of data points 
unchanged. In order to allow for non-circular motions both 
random and systematic, we adopt a constant uncertainty of 
7kms _1 for «t e rm- Hence each data point adds to xf n an 
amount 



^tcrm, model (0 ^term,data(0 

7kms~ 1 



(21) 



where w = 0.844 or w — 1.23 depending on whether / is 
greater than or less than 0. 

The rotation-curve data for R > Rq cannot be treated 
in an exactly analogous way because now two numbers con- 
tain significant uncertainties: d and W . Following Fich, Blitz 
& Stark (1989) we take the contribution to xt^t fr° m the *th 
data point to be 



w mm 

R>R 



In di — In d„ 



i(R) 



A In di 



+ 



W t -W n 



i(R) 



,(22) 
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Figure 2. Rotation curves of models 1-4 (thick solid lines). The 
circular velocities due the disk (dotted), bulge (short dashed), 
disk and bulge (dash dotted), and halo (long dashed) are also 
shown. 



where A hid and AW are the uncertainties in the observ- 
ables. In order to give each catalog the same weight, the data 
taken from Brand & Blitz (1993) and Pont et al. (1997) are 
weighted by w = 0.665 and 2.016, respectively, leaving the 
effective total number of data points unchanged at 119. To 
account for non-circular motions, a dispersion velocity of 
7kms _1 is quadratically added to the measurement errors 

Of Vlsr- 



5 RESULTS 

There are two aspects of any given model to consider: (i) how 
well does it fit the observational constraints, and (ii) how 
is its mass distributed. Since the observational constraints 
mostly relate to motions in the plane, these two questions 
are in large degree independent of one another. Table ^ lists 
the defining characteristics of each model and the quality of 
the fit to the data that it furnishes. Table ^ lists the values 
of each model's parameters. 

The most important parameter of the models proves to 
be Rd,*/Ro- The first four rows of Tables | and | refer to 
our 'standard models', which all have Rq — 8 kpc but adopt 
four values of R d .,/Ro. 0.25, 0.3, 0.35 and 0.4. Each model 
is determined by specifying the value of Rd,*/Ro and then 
solving for the values of the other model parameters which 
minimize_Yj ot . These values are given in Table ^| 

Fig. □ shows the circular-speed curves predicted by the 




Figure 3. Terminal velocities at R < Ro for models 1-4 and 2h. 
The data are from Weaver & Williams (1973, filled triangles), 
Bania & Lockman (1984, filled squares), Knapp et ai.(1985, open 
circles), and Kerr et ai.(1986, filled circles). 




R R„ 



Figure 4. W (equation |llj) versus R/Ro for models 1-4 (hardly 
to distinguish) and the two most extreme models 2b, i. The data 
are from Brand & Blitz (1993; triangles) and Pont et al. (1997; 
dots). 



four standard models together with the contributions to v c 
from each component. In all four models the circular speed 
declines outside Ro, although in Model 4 the decline is ex- 
tremely slow near Rq. 

As Rd,*/Ro increases from 0.25 to 0.4, the peak in 
the disk's contribution to v c moves outwards from 5 kpc 
and the amplitude of the disk's contribution to v c declines 
markedly. This decline in the disk's contribution to the in- 
ner circular-speed curve is compensated by an increase in 
the halo's contribution. This increase is achieved by making 
the halo more centrally concentrated, with the result that for 
Rd.*/Ro > 0.35 the halo does not differ greatly from a pure 
power-law component whose conbtribution to v c is nearly 
independent of radius. For Rd,*/Ro = 0.4 this contribution 
dominates v c at all radii. 

The parameter Ed, tot is the sum of the values for the 
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Table 3. Fit results 



Model 


«d,» 

n () 


Deviations from standard 


So 


A 


B 


Vc(Ro) 


2irG 


A^R<100kpc 


B W 


*?„ 


x 2 


X 2 th 




1 


0.25 


- 


43.3 


14.4 


-13.3 


222 


68.0 


6.56X10 11 


100 


91.2 


174 


1.92 


13.9 


2 


0.30 


- 


52.1 


14.3 


-12.9 


217 


72.2 


6.52x10" 


108 


82.3 


171 


0.53 


12.2 


3 


0.35 


- 


52.7 


14.1 


-13.1 


217 


72.5 


6.51X10 11 


111 


86.7 


175 


0.46 


12.5 


1 


0.40 


- 


50.7 


13.8 


-13.6 


220 


72.1 


6.04X10 11 


111 


93.7 


180 


0.94 


13.5 


2a 


0.30 


i?o = 7.5 kpc 


52.5 


15.1 


-12.2 


204 


71.3 


6.95x10" 


106 


83.5 


164 


0.70 


12.1 


2b 


0.30 


Ro = 8.5 kpc 


50.3 


13.8 


-13.3 


231 


71.9 


6.33X10 11 


108 


81.0 


176 


0.87 


12.7 


2c 


0.30 


7h = 1 


52.3 


14.3 


-12.7 


216 


72.0 


6.77X10 11 


105 


83.9 


171 


0.69 


12.4 


2d 


0.30 


7h = 1 and /3 h = 3 


53.3 


14.2 


-12.9 


217 


72.9 


6.33X10 11 


108 


84.3 


176 


0.66 


12.7 


2e 


0.30 


0.1 cos R/Rd added to exponent in (1) 


49.6 


14.3 


-13.0 


218 


71.3 


6.49x10" 


105 


85.5 


176 


0.78 


12.8 


2f 


0.30 


—0.1 cos R/Rd added to exponent in (1) 


52.8 


14.4 


-12.6 


217 


72.0 


6.77X10 11 


108 


80.6 


166 


0.42 


11.8 


2g 


0.30 


constraint agy/ = 140 ± 15kms _1 


51.4 


14.1 


-13.2 


218 


72.0 


6.39X10 11 


115 


95.1 


176 


2.99 


15.1 


2h 


0.30 


M h > 1.5 x 10 10 M Q 


44.1 


13.9 


-14.0 


223 


68.5 


6.01X10 11 


124 


131 


181 


1.71 


15.9 


2i 


0.30 


halo axis ratio q^ = 0.3 


40.0 


14.0 


-11.9 


207 


77.9 


5.59X10 11 


109 


88.8 


187 


2.68 


15.1 


4a 


0.40 


Ro = 7.5 kpc 


49.6 


14.3 


-13.6 


210 


71.5 


5.43x10" 


116 


103 


175 


1.16 


13.9 


4b 


0.40 


Ro = 8.5 kpc 


50.8 


13.4 


-13.7 


231 


72.4 


6.29x10" 


107 


89.1 


186 


1.50 


14.1 


4c 


0.40 


7h = 1 


49.8 


13.9 


-13.7 


221 


71.9 


5.75x10" 


106 


92.1 


178 


1.48 


13.8 


4d 


0.40 


7h = 1 and = 3 


48.3 


14.0 


-13.8 


222 


71.4 


3.73x10" 


107 


89.6 


173 


2.99 


14.7 


4e 


0.40 


0.1 cos ij/il d added to exponent in (1) 


51.7 


14.0 


-13.3 


218 


72.5 


6.36x10" 


111 


88.6 


179 


0.57 


12.9 


4d 


0.40 


—0.1 cos R/Rd added to exponent in (1) 


49.5 


13.7 


-13.9 


221 


71.7 


5.76x10" 


110 


99.1 


182 


1.41 


14.2 


lg 


0.40 


constraint itbw = 140 ± 15kms _1 


52.0 


13.8 


-13.4 


218 


72.4 


6.27x10" 


117 


106 


182 


2.87 


15.8 


4h 


0.40 


M b > 1.5 x 10 10 M Q 


48.3 


13.7 


-14.2 


223 


71.5 


4.89x10" 


120 


114 


181 


2.06 


15.5 


4i 


0.40 


halo axis ratio q^ = 0.3 


40.0 


13.6 


-12.3 


207 


80.2 


3.41x10" 


112 


105 


197 


5.57 


18.7 



The standard models 1 to 4 are determined by the choice of Rd */Ro (column 3), where Ro = 8 kpc. The non-standa rd m odels 2a-2i 
and 4a-4i are specified in column 4. Columns 5 to 11 give the best-fit values of the observables discussed in Section |3.3| . The values 
obtained for the x 2 s defined in Section ^ are given in the last four columns. The units are as usual: Ro in kpc; So, K Zy i,i/(2irG) in 
Mqpc -2 ; A, B in kms^kpc -1 ; -Mfi<iookpc hi Mq; and v c (Ro) and a Bw in kms -1 . 



three disks of the parameter Ed that is defined by equation 
(^). Table ^ shows that Ed.tot decreases by a factor of almost 
4 as Rd,*/Ro increases from 0.25 to 0.4. Consequently, the 
disk's contribution to v c at small radii decreases by a factor 
of 2, which is first compensated by a near doubling in the 
bulge mass, and then by a dramatic increase in the central 
density of the halo. Specifically, whereas for the smallest two 
values of Rd,,/Ro the halo has a hole at its centre (71, = —2), 
for all larger values of Rd t */Ro the halo density decreases 
outwards as ~ r~ 1,7 . 

The amplitude of the bulge's contribution to u c near 
the centre is largest for Rd,*/ Ro = 0.3 - smaller and larger 
values of this parameter yield bulges that are less massive 
by a factor in excess of 40 per cent. In all four models the 
velocity dispersion in Baade's window lies below the target 
value. As we indicated above, this shortfall probably reflects 
the fact that our models take no account of the elongation 
of the bulge along the line of sight. 

The model with the smallest value of Rd,*/Ro is nearly a 
maximum-disk model in the sense that there is a wide range 
of radii R 10 kpc within which nearly all of v c derives from 
the disk alone. In Model 2, which has Rd,*/Ro = 0.3, the 
disk and bulge together dominate v c at R ^ 10 kpc; the prin- 
cipal difference between Models 1 and 2 is the dominance of 
the bulge at R £ 2 kpc discussed above. For Rd^/Ro ii 0.35 
the disk's constribution to v c is nowhere as important as the 
halo, even though the disks of these models are only slightly 
less massive than those of models with Rd,*/Ro 0.3. 

All the standard models have masses Af^iookpc that 
lie below the target value. This is because the rotation-curve 
data for R > Ro force v c to decline just beyond the Sun. In 



extreme cases - for example Model 3 and Model 2a below 
- the constraint on -Mfj<iook P c obliges the halo to make an 
outwards-increasing contribution to v c even at the largest 
radii, with the result that v c is increasing at R — 100 kpc. 

All four standard models provide satisfactory fits to the 
constraints of Section 3.3 - see Table [| Fig. ^| shows that 
these models also provide excellent fits to the observed tan- 
gent velocities at 7? < Ro: the deviations between observed 
and predicted points are of the order of those expected to 
arise from spiral structure. Fig. ^ shows that all four stan- 
dard models provide essentially the same reasonable fit to 
the (widely scattered) rotation-curve data for R > Ro. 

Tables ^ and |i| also describe variants of the standard 
models. Models 2a and 2b explore the effect of changing 
Ro: reducing Ro from 8 kpc to 7.5 kpc increases A and re- 
duces Vc(Ro)- The top two panels of Fig. |^ show that it also 
changes the form of v c (R) from steady decline at all R > Ro 
to a steady rise at R ij 40 kpc. This coupling between the 
value of Ro and the rotation curve at R 3> Ro arises because 
the tangent velocities strongly constrain the mass distribu- 
tion at R < Rq. Hence any change in one component at 
R < Ro must be compensated by a change in another com- 
ponent. Since our model for the halo has only a few free 
parameters, a change in its density at R < Ro is accompa- 
nied by significant changes in density at R ^> Ro, and vice 
versa. 

Model 2c shows the effect of fixing the inner slope of 
the halo's density profile at 7h = 1; this value is motivated 
by Navarro, Frenk & White (1996), who find that in simula- 
tions of cosmological clustering, dark-matter halos tend to 
a universal profile that has a slope ~ 1 at small radii. Even 



© 1997 RAS, MNRAS 000, 000-000 



8 W. Dehnen and J. Binney 



Table 4. Best-fit values of the model parameters 



Model E d ,tot p ,t p ,h 7h A, r o,h n,h M d M h M h]<10kpc M hi<wokpc 



1 


1 QO^ 


n 4971 


7110 
U. t -L-1U 


2 


2 


. yjy 


o.oo 




^13 


U.O-lo 


2 81 


fiO 


2 


1 9ns 


n 7^fii 

U. / JU1 


1 9fi^ 


2 


2 


907 


1 flQ 
i.uy 




A 88 


Q1 7 


9 8Q 


A 

jy.^ 


Q 
O 


778 A 


u . o 


1 1 7Q 


1 8 
1 . o 


2 


009 


2 29 




4 46 






fiO ,; i 


■I 


536 


0.3 


0.2659 


1.629 


2 


.167 


1.899 


oo 


4.16 


0.364 


5.23 


55.9 


2a 


1222 


n 30^8 
u.oyoo 


u.uuuyiuu 


"1 8 
1 .O 


i 


-UD4: 


99 7Q 
zz. ( y 




4 32 




2 46 


(\A 7 


91 -> 
Z u 


1 lOO 


n 78 7A 


1 A9^ 


9 


9 
Z 


^ r i 3 


1 733 


OG 


^ 39 


u.yoo 


i 


^7 1 


2 c 


1213 




n 030^ 
u.uouoo 


1 


2 


. zuy 


U. lOO 




1 


fl 797 
U. 1 Z I 


9 QQ 
z.yy 


62 1 


2d 


1 938 






1 




Q 
O 


21 8 




K 
■ ) 


n 8i q 


9 83 
z.oo 


^7 ^ 


9 C 


1 9fifl 






2 


2 


425 


1 8fl9 




^ OA 


fl 7Q8 
u. ( yo 


9 88 
z.oo 




9f 
Zl 


± ±Z ± 


U.o 


D A1 7Q 


1 H 
1 .0 


i 


888 


1 


■OO 


1.0 


n 3ftzi 


3 7 
a. / 


ft9 7 
DZ. ( 


9o- 


1 1 Q3 


n q^3 


n 839fi 


.2 


2 


447 


1 Q1 Q 
i . y iy 




4 82 


11ft 


9 77 
Z. ( / 


^7 Q 

oi.y 


91 1 
Zil 


1 09A 


1 937 
1 . Zo i 


1 fl71 
l.U l 1 


Z 


9 
Z. 


888 

.000 


9 877 
Z.o 1 1 


OG 


A 1 3 


"1 r ! 

1 ,o 


O.'l 




2i 


928.2 


0.3380 


0.8514 


1.8 


i 


.868 


1 


OO 


3.75 


0.41 


2.92 


51.7 


4a 


526.5 


0.3 


0.7555 


1.757 


2 


.108 


1 


oo 


3.59 


0.364 


4.43 


50.4 


4b 


534.9 


0.3 


0.05205 


1.556 


2 


.386 


5.239 


oo 


4.7 


0.364 


6.26 


57.8 


4c 


527 


0.4147 


1.293 


1 


2 


.233 


1 


oo 


4.09 


0.503 


5.21 


52.9 


4d 


511 


0.6507 


0.1101 


1 




3 


5.236 


oo 


3.97 


0.789 


5.13 


32.5 


4e 


584 


0.3 


0.1078 


1.764 


2 


.076 


2.628 


oo 


4.45 


0.364 


4.74 


58.8 


4f 


491.3 


0.3 


0.3926 


1.5 


2 


-256 


1.764 


oo 


3.91 


0.364 


5.66 


53.4 


4g 


549.6 


0.3 


0.7447 


1.8 


2 


.066 


1 


oo 


4.27 


0.364 


5.04 


58.1 


Hi 


510.6 


1.237 


3.313 


-2 


2 


.672 


1.262 


oo 


3.97 


1.5 


4.44 


43.4 


4i 


423 


0.3 


2.799 


1.282 


2 


.336 


1 


oo 


3.29 


0.364 


3.95 


30.5 



Columns 2 to 8 give the values of the parameters obtained from the fitting procedure; surface and volume densities are given in 
Mq pc — 2 and Mq pc — 3 respectively. The last four columns give, in units of 10 10 Mq, the total masses of disk and bulge, and the 
mass of the halo within lOkpc and 100 kpc. 



though Model 2 has 7h = —2, fixing 7h in this way has a 
negligible effect on the observable properties of the model, 
and causes a negligible increase in Xt a t ■ This is because the 
halo makes a negligible contribution to the central density 
even for -yh = 1. Moreover, the increase in 7h is at small radii 
largely compensated for by a six-fold increase in the halo's 
break radius ro,h- 

Model 2d explores the effect of fixing both the inner and 
the outer slopes of the halo to values, 7h = 1,(3^ = 3, that 
are suggested by Navarro et al. (1996). Again the model's 
observables and value of change remarkably little, while 
ro,h increases again, this time by a further factor of ~ 3. 

Models 2e and 2f show the effect of employing a slightly 
non-exponential disk: a term ± 0.1 cos(i?//Zd) is added to 
the exponent in equation (Q) . The lower two panels of Fig. ^ 
show that at R ^ 10 kpc this extra term changes the balance 
between the contributions to ?; c of the bulge, disk and halo 
without significantly affecting the overall rotation curve. At 
large R the effect of the additional term is more dramatic in 
that with a plus sign v c declines steadily at R > Ro, while 
with a minus sign it rises gently at the largest radii. Again 
we encounter the consequence of the rotation curve at large 
radii being inadequately constrained by the observations: it 
mirrors changes in the functional form of the halo that are 
designed to fit the data at R ^ Ro. 

Model 2g shows that increasing <tbw, the target dis- 
persion in Baade's window, by 23kms _1 has little effects 
on the derived model: the predicted dispersion is grudgingly 
raised from 108kms _1 to 115kms _1 , and Xtot increases sig- 
nificantly. 

Microlensing surveys have suggested that the bulge may 
be significantly more massive than studies of the tangent 
velocities would imply (e.g., Bissantz et al. 1996). Model 



Table 5. Parameter Correlations for models 2 (upper right) and 
4 (lower left) 





P0,h 


7h 


/9b 




n,h 


Po,b 


Ed, tot 


P0,h 




-0.999 


-0.953 


-0.999 


-0.909 


0.917 


-0.161 


7h 


-0.714 




0.952 


0.998 


0.907 


-0.920 


0.153 




-0.990 


0.646 




0.968 


0.989 


-0.788 


0.112 




-0.999 


0.705 


0.992 




0.929 


-0.900 


0.153 


r t ,h 


-0.963 


0.643 


0.984 


0.965 




-0.729 


0.109 


P0,b 


-0.877 


0.312 


0.891 


0.882 


0.846 




-0.175 


^d, tot 


0.089 


-0.126 


-0.061 


-0.089 


-0.005 


-0.013 





2h shows the effect of imposing the constraint Mb u i gc > 
1.5 x 10 10 M Q . This causes v c to be dominated by the bulge 
out to ~ 3.5 kpc. It also, by the mechanism described above, 
changes the form of v c (R) at large R so that with a more 
massive bulge v c is predicted to decline significantly more 
steeply at 7? ^ 50 kpc. Like the other model with an arti- 
ficially enhanced bulge (Model 2g), this model has an ana- 
malously large value of Xjf ot . 

It is not clear that a massive halo should be nearly 
round (axis ratio qh = 0.8) as in all the models described 
above, so Model 2i explores the effect of assuming that the 
halo is highly flattened: axis ratio — 0.3. This has the 
effect of making the halo virtually a pure power-law compo- 
nent that is dynamicaly dominant at all radii. Consequently, 
it predicts that the rotation curve rises at the very largest 
radii. The massive, flattened halo contributes strongly to 
K Zl i.i, so that the target value of A" 2 ,i.i is significantly over- 
shot, despite the local disk surface density, Eo, being pushed 
right down to its floor value, Eo = 40 Mq pc -2 . For this mo- 
del Xtot ls on the high side. 
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Figure 5. Rotation curves of models 2a,b,e,f (thick solid lines). 
The circular velocities due the disk (dotted), bulge (short dashed), 
disk and bulge (dash dotted), and halo (long dashed) are also 
shown. 



Models 4a - 4i explore the effect on Model 4 of the 
changes that were made to Model 2 in making Models 2a - 
2i. Qualitatively the results are similar, but quantitatively 
they tend to be smaller because the halo is very much more 
important in Model 4 than it is in Model 2 and the effects 
of changes in the disk and bulge are relatively minor. 

Inspection of Table ^ shows that there is a general anti- 
correlation between the steepness of the halo's central den- 
sity profile and the bulge mass. In particular, when the halo 
density is strongly cusped at the centre (% > 1), the bulge 
mass is small: Mb ~ 0.37 x 1O 1O M0. Conversely, with the 
exception of the model with the smallest value of Rd^/Ro, 
the bulge mass is large (Mb 0.8 x 10 10 M ) when the halo 
has a hole at its centre (7h = —2). 

The asymptotic slope of the halo profile at large R is 
always larger than (3^ = 1.63 and no model has a finite halo 
cut-off radius r t ,^. That is, the data imply only that the halo 
ends at R > 100 kpc. 

The only models that comes near to violating the lower 
limit ( [rij ) on the local column density are Model 1, which 
has the smallest disc scale-length, Model 2h, which is obliged 
to have a massive bulge, and Models 2i and 4i, which have 
highly flattened halos. 

Table ^ gives the correlation matrix of the fitted param- 
eters for both Models 2 and 4: the upper right triangle is for 
Model 2 while the lower-left triangle is for Model 4. The 
parameters of the halo are strongly correlated with one an- 



other. Indeed, the halo's density normalization po,h, central 
slope 7h and scale length ro,h have almost unit correlations 
between them. The density normalization of the bulge, po^ 
is strongly correlated with all the halo parameters, espe- 
cially po,h, 7h and ro.h- By contrast, Ed, tot is only weakly 
correlated with the other parameters. These strong corre- 
lations between parameters reflect the existence of a long, 
level valley in parameter space. 



6 CONCLUSIONS 

We have fitted a multiparameter mass model to the avail- 
able kinematic data for the Milky Way. The wide variety of 
models that emerge from this fitting process demonstrates 
that the mass distribution within the Milky Way is currently 
ill determined. 

The principal problem in determining the Galaxy's 
mass distribution is that very few hard facts are available 
regarding the vertical distribution of the Galaxy's mass. 
Also the circular speed is observationally much less well con- 
strained at R > Ro than it is at R < Ro. These deficiencies 
in the spatial coverage of the data have several unfortunate 
effects. First they oblige use to represent the Galaxy as a su- 
perposition of components, and to adopt simple functional 
forms for the distribution of mass within each component. 

The assumed distribution of mass within the halo plays 
a particularly important and confusing role. Since we know 
nothing about the halo except what can be gleaned from 
studies such as this, we have allowed the halo density as 
much freedom as is compatible with (i) its being much 
thicker than the disk (axis ratio >= 0.3), (ii) its density 
function containing only a few free parameters, and (iii) its 
not being implausibly sharp-edged ( — 2 < 7h < 1, Ai > 1). 

With these assumptions we find that even the circular- 
speed curve of the best-fitting model depends significantly 
on both the adopted values of parameters such as Ra,*/Ro 
and Ro, and on the adopted functional forms of the compo- 
nents. The density model is even less well constrained by the 
data. In particular, remarkably small changes in v c (R) and 
the other observational constraints can be associated with 
dramatic changes in the distribution of mass between the 
different components, and the degree of central concentra- 
tion of the halo. 

A particularly disturbing phenomenon is that a change 
in a component that contributes to v c only at small radii 
causes the predicted value of v c to change at large radii. 
This connection between small and large scales within the 
Galaxy is established by the halo component, which must be 
allowed to contribute to the density at all radii and yet be de- 
termined by a small number of parameters. In general there 
is a clear need to model components non-parametrically or 
at least by functions that contain more parameters. How- 
ever, we do not yet have enough observational constraints 
to constrain adequately models that are significantly more 
complex than those used here. 

As is traditional, we have represented the Galaxy as 
a superposition of components that individually represent 
plausible stellar systems. The justification for the use of 
such components in preference to a family of orthogonal 
(and therefore non-positive) functions, is the hypothesis that 
these components do, in fact, represent real physical sys- 
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tems. While there can be no doubt of the reality of the 
disk and bulge, the status of the halo is entirely speculative. 
Indeed, although observations of external galaxies clearly 
require high mass-to-light ratios at large radii, it does not 
follow that these reflect the existence of a physically distinct 
dark halo; it is perfectly possible that the mass-to-light ra- 
tio of the disk or bulge increases strongly away from the 
Galactic centre. If the halo does represent a distinct dyna- 
mical entity, then in an exercise like the present one it should 
emerge with a dynamically plausible density profile. Table [l| 
shows that it frequently fails this test: in just over half the 
models, the central density slope 7h lies at one or other ex- 
treme of its permitted range —2 < 7h < 1.8. In four models 
(Models 3, 2a,f,i) the halo's density profile is an essentially 
featureless power law of slope ~ —1.8. In seven models the 
halo has a hole at its centre. JV-body simulations of structure 
formation offer no encouragement to the idea that the Ga- 
laxy's distribution of axions or other exotic particles would 
have a hole at its centre. Nor do they suggest that it should 
be a featureless power law (Navarro et al. 1996). 

The road to an improved understanding of the Galaxy's 
mass distribution must lie with the introduction of more ob- 
servational constraints, especially ones that relate to R > Rq 
and z <; R. The goal of later papers in this series is to bring 
to bear on this problem observations of halo stars, which 
should provide a wealth of information about the density at 
z <; R, and, to a lesser extent, about the density at R > Ro. 
The models described here simultaneously provide starting 
points for this enterprise and demonstrate its urgency by 
underlining that at the present time we know depressingly 
little about the distribution of mass within our own Galaxy. 
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